X = LoadData('data4.bin');
[X_norm,t] = ecg_normalize(X);
Samp_freq=553.625;
Rs=20;
Fp= [0.5 96];


        Wp = Fp / Samp_freq *2;
        Ws = [eps 150] / Samp_freq*2;
        Rp = 1;

        %generate cheby1 filter
        [n Wc] = cheb1ord(Wp,Ws,Rp,Rs);
        [b,a]=cheby2(n,Rp,Wc);

        %filter signal
        filtfilt(b,a,X_norm);
        


Depth=10;
Width=1;
F=50;

Ws = [F-Width/2 F+Width/2] / (Samp_freq/2); % Stop band
Wp = [Ws(1)*0.7 Ws(2)*1.3]; % Pass band
Rp=1;  %Maximium ripple in the pass band (notch) in dB
Rs = Depth; %Minimium attenuation in the stop band (notch) in dB

[n,Wc] = buttord(Wp,Ws,Rp,Rs);
[b1,a1] = butter(n,Wc,'stop');
filtfilt(b1,a1,X_norm);

Ws2 = [100-Width/2 100+Width/2] / (Samp_freq/2); % Stop band
Wp2 = [Ws2(1)*0.7 Ws2(2)*1.3]; % Pass band
Rp2=1;  %Maximium ripple in the pass band (notch) in dB
Rs2 = Depth; %Minimium attenuation in the stop band (notch) in dB

[n,Wc2] = buttord(Wp2,Ws2,Rp2,Rs2);
[b2,a2] = butter(n,Wc2,'stop');
filtfilt(b2,a2,X_norm);

a3 = conv(a, a1);
b3 = conv(b, b1);
a4 = conv(a3, a2);
b4 = conv(b3, b2);

fvtool(b4,a4);